RDは計量経済学で用いられる手法の一つです。
人道的、社会的制約のためにランダム化比較実験が行えない状況を想定します。
そのような状況で因果効果を算出するために、処置の割り当てが擬似的にランダムとなる状況を利用する一連の手法が存在します。
RDはその代表的な手法であり、ImbensやAngristの手により近年急速に発展しました。
この記事では、RDの基本的な発想を解説しつつ、Juliaを使って実際にRDによる因果効果の測定を行ってみます。
RDの理論的側面を解説します。
RDには_Sharp RD_と_Fazzy RD_の二種類が存在します。ここではSharp RDのみ解説します。（Fuzzyはそれがわかれば簡単なので、教科書を読んでください） 特に、angristなどでさっと流されるband width selectionについてそのやり方を書いたので読んどいてください。
次のように具体的な状況を考えます。
経済学者の河原崎は東大に行くことが学力向上に寄与するかどうかに関心を持っています。 彼は学力の客観的な指標の作成に成功しており、そのテストをPikaと名付けました。Pikaは特殊なテストで、得点分布は正規分布に従うことがわかっています。
この夏、河原崎は東大とW大学の学生に対して大規模な学力調査を行い、各大学から大学4年生1000人分のPika点数データを入手しました。
そのデータは以下のような構造をしていました。
In [2]:
# パッケージインストール
using Gadfly
using DataFrames
using GLM
In [3]:
srand(22)
n = 1000
todai = Array{Float64}([i + 1 for i in randn(n)])
wuni = Array{Float64}([i - 0.5 for i in randn(n)])
plot(
layer(x = todai, Geom.histogram),
layer(x = wuni, Geom.histogram, Theme(default_color=colorant"red")),
Guide.xlabel("Pika Score"), Guide.ylabel("Population"), Guide.title("Pika Score")
)
Out[3]:
赤がw大学のpika scoreで、青が東大のpika scoreです。
河原崎はこのデータを見て、「東大に行くと学力が向上しそうだ！」と思いました。
彼はこの仮説を検証するために、学生時代に習った回帰分析を使ってpika scoreを東大生のダミー変数と定数項に回帰してみることにしました。
回帰式は下のようになります。
結果は以下のようでした。
In [4]:
srand(22)
dummy = ones(2n)
dummy[1001:end] = 0
pika = [todai; wuni]
X = reshape([ones(2n) dummy], 2n, 2)
slope_coef = (inv(X' * X) * X') * pika
coef = ("alpha", "beta")
for i in 1:length(coef)
c = coef[i]
d = slope_coef[i]
println("$c = $d")
end
t値等も合わせて結果を表記すると以下のようになります。
「東大ダミーは有位に正だ！やっぱり東大に行くと学力が向上するんだ！」と河原崎は思いました。
しかし、経営学者の正田がこれにまったをかけます。
「東大生はもとからW大学の学生よりも頭がいいんだから、その結果は東大に入ったことによる学力向上が測定されているとは言えないんじゃない？」
たしかにその通りです。 しかし、河原崎はこれ以上の良い方法を知りません。
どのようにすれば彼らは東大に入ることの効果を識別することができるでしょう？
2013年度東京大学文科Ⅲ類の合格最低点は347.2111点でした。つまり、347.211点を0.001点でも上回った受験者は東大に合格し、0.001点でも下回った受験者は東大には行くことはできませんでした。
RDはこの事実を使って、347.211点をわずかに超えて合格した人と、わずかに下回って不合格になった人とで従属変数を比較する手法です。
今、A君の合計得点が348点で、B君の合計得点が347点だとしましょう。
A君はめでたく東大に合格できましたが、B君は残念ながら不合格、W大学に通うことになりました。
さて、この二人について先ほど正田が述べた理屈は当てはまるでしょうか？
二人の得点差はわずかに1点。これはセンター試験の国語で1問落としたとか、世界史でたまたま直前に見てたところが出たとか、そのぐらいの偶然で十分に生まれる差です。
しかし、このわずか1点によって彼らが後の4年間を過ごす大学が変わってくるわけです。
ここで東大に入ることの学力向上への効果を見たいわけですから、東大入学を処置と考えると、 ほとんど能力に差のない二人に処置の割り当てがランダムに行われる、という状況が存在していることがわかります。
これがまさに自然実験であり、擬似的なランダム化比較実験が行われていると解釈できる状況です。
Angristでやった通り、ランダム化が行われていれば処置は潜在変数と独立となるため、処置群と対照群とで従属変数の平均値を比較すれば、それを平均的な因果効果として算出することができます。
今回の文脈で言えば、A君とB君のPika Scoreの点を比べれば、それは東大で勉強したことによる因果効果を算出している、ということです。
以下で実際のデータを模して作った擬似データを使ってこの状況をグラフにします。
In [5]:
# 前期試験の得点分布作成
srand(22)
cut = 347.211
n_1 = 500
n_2 = 300
n_3 = 150
n_4 = 50
n_array = [n_1, n_2, n_3, n_4]
pop_todai = Array(Float64)
for (i, n) in enumerate(n_array)
pop_todai = [pop_todai ; [j*25 + cut + 25*(i-1) for j in rand(n)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])
pop_w = Array(Float64)
for (i, n) in enumerate(n_array)
pop_w = [pop_w ; [cut - j*25 - 25*(i-1) for j in rand(n)]]
end
pop_w = Array{Float64}(pop_w[2:end])
plot(
layer(x = pop_todai, Geom.histogram),
layer(x = pop_w, Geom.histogram, Theme(default_color=colorant"red")),
Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Population")
)
Out[5]:
In [6]:
srand(22)
# 前期試験の得点とpika scoreの散布図作成
s_1 = reshape(pop_todai,1000,1)
s_2 = reshape(pop_w,1000,1)
score = vcat(s_1, s_2)
data = hcat(X, pika)
data = hcat(data, score)
n = 1000
cutoff = 347.211
# 回帰直線
d_todai = DataFrame(X = pop_todai, Y = todai)
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = pop_w, Y = wuni)
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)
# plot
todai_layer = layer(x = pop_todai, y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = pop_w, y = pre_w, Geom.line, Theme(default_color = colorant"black"))
gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of attend to todai on pika score is $gap")
plot(x = score, y = pika, todai_layer, w_layer, xintercept=[347.211],
Scale.x_continuous(minvalue= 200, maxvalue=500),
Geom.vline(color = colorant"orange"),
Geom.point,
Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
Theme(default_point_size = 1.5pt)
)
Out[6]:
上の散布図をみてください。
オレンジの線が合格最低点である347.211点です。
オレンジの線を境にして、左側では低めに、右側では高めにpika scoreが出ていることが見て取れるでしょう。先ほど説明した状況になっているのがわかりますか？
ここでは東大に入学するという処置を受けられるか否かの境でpika scoreが不連続となっているのですね。
後の手順は簡単です。
オレンジ線の左側のデータのみで線形回帰を行いcut off(ここでは347.211点)における予測値を算出し、右側においても同様の手順でcut off における予測値を算出します。その時の予測値の差が求めている因果効果です。
ただ、一点注意しなくてはならない点があります。
それはデータをどこまで使うのか、という問題です。
知りたい因果効果が現出するのがオレンジ線上であることを考えると、前期試験の点数がとても高い学生や、極端に低い学生のデータは今必要ではありません。
私たちが知りたいのは、左側で作った予測値のオレンジ線における極限値と右側で作った予測値のオレンジ線における極限値の差です。
それゆえにオレンジ線を含む短い区間のデータのみに注目して切片を予測するのは理にかなったことといえます。
しかし、ある程度遠くのデータまで含めないと推定値が安定しないこと、またあまりにlocalな結果となってしまうことも事実です。
RDにおいては、このジレンマの中で最適なband width、すなわちデータを用いる幅を決定することが肝となってきます。
その際に理論的な背景を与えるのがLocal Linear Regressionです。以下ではこのLocal Linear Regressionについて解説し、上のデータを少し変えて最適なbandwidthを実際に求めてみたいとおもいます。
上のデータは面白くないので、相関をつけます。さらに人数を実際に近づけて合格者400人、不合格者400人にします。 以下のモデルでpika scoreを算出し、上と同様のプロットを行います
In [7]:
srand(22)
# データ作成
srand(22)
n = 400
cutoff = 347.211
n_1 = n/2
n_2 = 3n/10
n_3 = 3n/20
n_4 = n/20
n_array = Int64[n_1, n_2, n_3, n_4]
pop_todai = Array(Float64)
for (i, r) in enumerate(n_array)
pop_todai = [pop_todai ; [j*25 + cutoff + 25*(i-1) for j in rand(r)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])
pop_w = Array(Float64)
for (i, r) in enumerate(n_array)
pop_w = [pop_w ; [cutoff - j*25 - 25*(i-1) for j in rand(r)]]
end
pop_w = Array{Float64}(pop_w[2:end])
s_1 = reshape(pop_todai,n,1)
s_2 = reshape(pop_w,n,1)
score = vcat(s_1, s_2)
newpika = Array(Float64, 2n)
for i in 1:2*n
if i <= n
newpika[i] = 10*log(score[i]) + randn()
else newpika[i] = 9.8*log(score[i])+ randn()
end
end
# 回帰直線
d_todai = DataFrame(X = score[1:n], Y = newpika[1:n])
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = score[n + 1:end], Y = newpika[n + 1 : end])
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)
gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of the attend to todai on pika score is $gap")
# plot
todai_layer = layer(x = pop_todai, y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = pop_w, y = pre_w, Geom.line, Theme(default_color = colorant"black"))
plot(x = score, y = newpika, todai_layer, w_layer, xintercept=[347.211],
Scale.x_continuous(minvalue= 200, maxvalue=500),
Geom.vline(color = colorant"orange"),
Geom.point,
Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
Theme(default_point_size = 1.5pt)
)
Out[7]:
local linear regressionは、その名の通り一部のデータを使って線形回帰を行う手法です。RDの文脈では、cut off周囲のデータのみを使って左側と右側で二つの線形回帰を行います。推定された係数をもとにcut offにおける推定値の差を算出し、それを持って処置の効果とする手法です。
その際にどこまでを周囲とするか、すなわちband widthを決めなくてはなりません。一般にこれは恣意的に決めるものではなく、data drivenな手法で求めるものされています。
band widthの決定には主に二つの手法が存在します。
1つ目はRule of Thumbで初期のband widthを決定し、それを改善していくもの。
2つ目はcross validationを行って、平均二乗誤差を最小にするband widthを選択するものです。
今回は2つ目の手法を紹介します。
cross validationとは一般に以下のような手法を指します。
Local Linear Regressionの文脈では、これは以下の用に用いられます。
では、以下で実際にband widthを選択してみましょう。
In [65]:
# データ作成
srand(22)
n = 400
cutoff = 347.211
n_1 = n/2
n_2 = 3n/10
n_3 = 3n/20
n_4 = n/20
n_array = Int64[n_1, n_2, n_3, n_4]
pop_todai = Array(Float64)
for (i, r) in enumerate(n_array)
pop_todai = [pop_todai ; [j*25 + cutoff + 25*(i-1) for j in rand(r)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])
pop_w = Array(Float64)
for (i, r) in enumerate(n_array)
pop_w = [pop_w ; [cutoff - j*25 - 25*(i-1) for j in rand(r)]]
end
pop_w = Array{Float64}(pop_w[2:end])
s_1 = reshape(pop_todai,n,1)
s_2 = reshape(pop_w,n,1)
score = vcat(s_1, s_2)
newpika = Array(Float64, 2n)
for i in 1:2*n
if i <= n
newpika[i] = 10*log(score[i]) + randn()
else newpika[i] = 9.8*log(score[i])+ randn()
end
end
# cross validation
band_widthes = collect(1.0:0.5:15.0)
validation_l = Array{Float64}(length(band_widthes))
validation_r = Array{Float64}(length(band_widthes))
data = DataFrame(Y = newpika, X = squeeze(score',tuple(1)))
for (t, i) in enumerate(band_widthes)
gaps = Float64[]
temp_d = data[cutoff .> data[2] .> cutoff-i, :]
temp_f = data[cutoff .> data[2] .> cutoff-i, :]
for j in 1:size(temp_d)[1]
temp_d[j, 1] = NA
independent = temp_f[j, 2]
dependent = temp_f[j, 1]
temp_res = lm(Y ~ X, temp_d)
predictval = predict(temp_res, [1, independent]')
squared_gap = ((dependent - predictval)[1])^2
push!(gaps, squared_gap)
temp_d[j, 1] = dependent
end
validation_l[t] = mean(gaps)
end
for (t, i) in enumerate(band_widthes)
gaps = Float64[]
temp_d = data[cutoff+i .> data[2] .> cutoff, :]
temp_f = data[cutoff+i .> data[2] .> cutoff, :]
for j in 1:size(temp_d)[1]
temp_d[j, 1] = NA
independent = temp_f[j, 2]
dependent = temp_f[j, 1]
temp_res = lm(Y ~ X, temp_d)
predictval = predict(temp_res, [1, independent]')
squared_gap = ((dependent - predictval)[1])^2
push!(gaps, squared_gap)
temp_d[j, 1] = dependent
end
validation_r[t] = mean(gaps)
end
leftband = band_widthes[indmin(validation_l)]
rightband = band_widthes[indmin(validation_r)]
for (i, n) in enumerate(band_widthes)
println("mean squared error for left when band = $n", " is ", validation_l[i])
println("mean squared error for right when band = $n", " is ", validation_r[i])
end
println("selected band width for the left side is $leftband")
println("selected band width for the right side is $rightband")
上の例では cutoffの左側ではbandwidth = 2.5が選択され、右側ではbandwidth = 5.0が選択されました。ちなみにNaNは、その範囲にデータが存在しなかったことを意味しています。
では、上記の最適なbandwidthに従ってデータを以下にプロットし、先と同様に東大進学の効果を計算してみましょう。
In [63]:
srand(22)
# 回帰直線
cutoff = 347.211
n = 400
d_todai = DataFrame(X = score[1:n], Y = newpika[1:n])
d_todai = d_todai[cutoff + rightband .>d_todai[1], :]
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = score[n + 1:end], Y = newpika[n + 1 : end])
d_w = d_w[d_w[1] .> cutoff-leftband , :]
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)
data = DataFrame(Y = newpika, X = squeeze(score',tuple(1)))
gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of the attend to todai on pika score is $gap")
# plot
todai_layer = layer(x = d_todai[1], y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = d_w[1], y = pre_w, Geom.line, Theme(default_color = colorant"black"))
plot(x = data[cutoff+rightband .>data[2].> cutoff-leftband, 2], y = data[cutoff+rightband .>data[2].> cutoff-leftband, 1], todai_layer, w_layer, xintercept=[347.211],
Scale.x_continuous(minvalue= 340, maxvalue=350),
Geom.vline(color = colorant"orange"),
Geom.point,
Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
Theme(default_point_size = 1.5pt)
)
Out[63]:
800人全てのデータを使って算出した時よりも効果は大きく推定されました。
In [ ]: